% EDITED BY FRANCESCO MARANGONI ON 2/8/2011
% CALCULATES INST. VELOCITY AND MSD CONSIDERING DISTANT JOINING
% CALCULATES ARREST COEFFICIENT
% DOES NOT HANDLE SIGNALING

function [varargout] = TrackerRead_March(infile, cyc_time, max_vel, min_track_length, min_disp, min_inst_vel, varargin)
% definisce il file come file di funzione. Variabili output tra[], variabili input tra()
% infile e' il file output di Imaris, cyc_time, max_vel, min_track_length, min_disp, ed altre variabili (varargin) provengono dall'interfaccia grafica di CellTracker

wh = waitbar(0, 'Reading and Analyzing file...');

fid = fopen(infile);
% apre la variabile "infile" assieme al file identifier (funzione fopen, -1 non si puo', 1 si puo', 2 standard error)
if (fid <= 0),
    error('Cannot open data file provided.  Check name and/or path.');
end  % Termina sempre un ciclo

all_names = {};
all_data = {};
data = [];
num_tracks = 1;
numlines = 1;
current_track = NaN;
%{} e' cella vuota di un array di array, [] e' array vuoto

%LEGGE FILE INPUT COME ARRAY DI ARRAY
while ~feof(fid),   % finche' (~) NON (feof) fine del file identifcato da fid 
    tline = fgetl(fid); %legge una linea senza il terminatore (fgetl) e la assegna a tline
    if (length(tline) < 4), % se la lunghezza degli elementi riga e' meno di 4
        %ignore this bad bad line
    elseif (strcmp(upper(tline(1:4)), 'NAME')),  % o se gli elementi 1:4 di tline sono NAME (uppercase)
        % We are in the first line.
    else  %in tutti gli altri casi
        lineofdata = {};
        datacount = 1;
        while (length(tline) > 0),  %finche' la linea letta (tline) ha elementi (e' >0)
            [token, tline] = strtok(tline, char(9));  %attribuisci a token il valore della linea fino a tab (char(9), e salva il resto in tline
            lineofdata{datacount} = token; %attribuisci all'array di array lineofdata nella posizione 1 (poi 2, 3, e cosi' via per ogni campo sono 6 nel file input) il valore di token (cioe' del campo stesso)
            datacount = datacount+1; %incrementa datacount (il valore di campo)
        end  %questo ciclo memorizza una linea di dati suddividendo ogni riga di dati in campi in base al carattere tab
        
        if (strcmp(lineofdata{2}, 'N/A')),  % se la pos 2 di lineofdata e' N/A
            %ignore
        elseif (str2double(lineofdata{2}) ~= current_track),  % o se la pos 2 di lineofdata (convertita in numero double precision) e' diversa (~=) dal valore di current_track
            % We are at the beginning of a cell track
            if ~isempty(data),  %se l'array data non e' vuoto
                all_data{num_tracks} = data;  %attribuisci all'array di array all_data (alla posizione num_tracks, inizialmente 1) il valore di data
                num_tracks = num_tracks + 1;  %incrementa num-tracks
            end
            %questo ciclo suddivide l'input in tracce, aggiungendo riga per riga il risultato all'array di array all_data

            current_track = str2double(lineofdata{2});  % stabilisce il numero di traccia corrente
            numlines = 1;  %resetta numlines a 1
            data = [];  %resetta data
            for i = 1:6,  %Modified 2/19/11 to handle cell legnth as input
                data(numlines, i) = str2num(lineofdata{i+1});  %crea il vettore "data" dal vettore complessivo "lineofdata" senza la posizione 1 (cioe' il nome dell'oggetto) e fa partire una conta progressiva (numlines) degli oggetti che compongono ogni traccia.
            end
            numlines = numlines + 1; %incrementa il numero di linea
            
            lineofdata = {};  %resetta lineofdata, diventato inutile perche' ha passato i dati (tranne il nome dell'oggetto) a "data"

        else
            %Data past first line.
            for i = 1:6,  %Modified 2/19/11 to handle cell legnth as input
                data(numlines, i) = str2num(lineofdata{i+1});
            end
            numlines = numlines + 1;
        end
	% stessa cosa per le righe successive. Finora, distingue le tracce solo sulla base del numero progressivo
    end
end
% Questo blocco legge il file input suddividendolo in campi (con l'eccezione del nome oggetto, che viene eliminato) e tracce!!

if ~isempty(data),
    all_data{num_tracks} = data;
    num_tracks = num_tracks + 1;
end
% salva l'ultima traccia

fclose(fid);    %chiude il file input
waitbar(0.25);  %fa procedere la waitbar ad un quarto

% open output files
analysis_fid = fopen(sprintf('%s_analysis.txt', strtok(infile, '.')), 'w');
if (analysis_fid <= 0),
    error('Cannot create analysis output file.');
end

summary_fid = fopen(sprintf('%s_summary.txt', strtok(infile, '.')), 'w');
if (summary_fid <= 0),
    error('Cannot create summary output file.');
end

msd_fid = fopen(sprintf('%s_MSD.txt', strtok(infile, '.')), 'w');
if (msd_fid <= 0),
    error('Cannot create MSD output file.');
end

%try
    % print header con tabulazioni (\t) e ritorni di linea (\n)
    fprintf(analysis_fid, 'Track #\tStep #\tCentroid X\tCentroid Y\tCentroid Z\tDelta 3D Dist (microns)\tDelta x\tDelta y\tDelta z\tInst. 3D Vel (microns/min)\tInst. Vel. X\tInst. Vel. Y\tInst. Vel. Z\t3D Disp (micron)\tDisp. X\tDisp. Y\tDisp. Z\tAngle Change\tInst. Vel. 2D\tMotile Angle Change\tMotile Angle Change 2\tTimepoint\tDelta Time\tCell length\n');
    fprintf(summary_fid,'Track #\tTotal Dist.\tTotal X\tTotal Y\tTotal Z\tMean Vel. 3D (microns/min)\tAvg. Vel X\tAvg. Vel Y\tAvg. Vel Z\tTotal Disp\tTotal Disp. X\tTotal Disp. Y\tTotal Disp. Z\tAvg. Angle Change\tConfinement Ratio\tAvg. Vel 2D\tTrack Duration (min)\tArrest Coefficient\tAvg. Cell length\tExceeded Max Allowed Velocity\n');
    fprintf(msd_fid,'Track #\tStep Size\tD(t)\tMSD\n');
    
    for i=1:length(all_data),  %i diventa il numero di traccia
        cur_data = all_data{i};
        
        if (size(cur_data, 1) < (min_track_length - 1)),  
            continue;
        end 
	%se la lunghezza della traccia (cur_data, campo 1) e' inferiore o uguale al minimo consentito, salta l'analisi      
        
        % distances
        del_pos = diff(cur_data(:, 3:5));  %attribuisce al vettore del_pos la differenza (diff) tra due posizioni consecutive (:) dell'array current data alle posizioni 3,4,5 (x,y,z)
        del_dist = sqrt(del_pos(:,1).^2 + del_pos(:,2).^2 + del_pos(:,3).^2); %attribuisce al vettore del_dist i valori di distanza.
        del_dist_2D = sqrt(del_pos(:,1).^2 + del_pos(:,2).^2); % line added 2/9/2011
	
	% computation of time distance between points ADDED 1/16/2011
	del_time = diff(cur_data(:, 2));   
        
        % velocities
        inst_vel = del_pos./(cyc_time/60);  % ORIGINAL LINE il vettore inst_vel riporta le velocita' istantanee su x,y,z SBAGLIA SE DUE PUNTI NON SONO CONSECUTIVI.
        inst_vel_3D = (del_dist./del_time)/(cyc_time/60);  % MODIFIED LINE 1/16/2011
        inst_vel_2D = (del_dist_2D./del_time)/(cyc_time/60);   % MODIFIED LINE  2/9/2011
        
        % displacements
        del_disp = cur_data(2:end, 3:5) - repmat(cur_data(1, 3:5), [size(cur_data, 1)-1 1]);
        del_disp_3D = sqrt(del_disp(:,1).^2+del_disp(:,2).^2 + del_disp(:,3).^2);

	% attempts to calculate a signaling index based on volume, sphericity, intensity  FRANCESCO MARANGONI 1/25/2011
	% sign_index = 1000.*(cur_data(:,7).*cur_data(:,8))./(cur_data(:,6).^2);  % signaling index = sphericity*intensity/volume
    % sign_index = cur_data(:,9);  %  here considers manual annotation of signaling
        
        % angle change
        del_angle = [];
        for j = 1:(size(del_pos,1)-1),
            if (norm(del_pos(j,:)) == 0) || (norm(del_pos(j+1,:)) == 0),
                del_angle(j) = 0;
            else
                del_angle(j) = acos(dot(del_pos(j,:), del_pos(j+1,:))/(norm(del_pos(j,:))*norm(del_pos(j+1,:))))*180/pi;
            end
            % A positive cross product along z is counter-clockwise
            %  We shall call this a negative angle
            temp = cross(del_pos(j,:), del_pos(j+1, :));
            del_angle(j) = -sign(temp(3))*del_angle(j);           
        end
   
        % make the column of angle changes without non-motile cell angles
        % (defined by comparing inst. disp. to user's entered min)
        motile = [];
        for c = 1:length(del_disp_3D)-1,
            if (del_disp_3D(c+1) >= min_disp),
                motile(c) = del_angle(c);
            else
                motile(c) = NaN;
            end
        end

	% ARREST COEFFICIENT COMPUTATION (CONSIDERS DISTANT JOINING)
	% FRANCESCO MARANGONI 1/17/2011
	nonmotile_arr_coeff = [];
	for d = 1:length(inst_vel_3D),
	   if (inst_vel_3D(d) <= min_inst_vel),
		 nonmotile_arr_coeff(d) = del_time(d);
           else
             nonmotile_arr_coeff(d) = 0;
           end
        end
        
        
        
        exceeded = find(abs(inst_vel_3D) >= max_vel);  %trova indici delle velocita' istantanee superiori al limite prefissato
        exceeded = [0; exceeded; size(inst_vel_3D,1)+1];  % FONTE DEL PROBLEMA: SEPARA TRACCE SULLA BASE DEL SUPERAMENTO DELLA VELOCITA' MAX CONSENTITA!!!!!!!!!
        
        for j = 1:(length(exceeded)-1),
            start_pos = exceeded(j) + 1;
            end_pos = exceeded(j+1) - 1;
            if ((end_pos - start_pos) < (min_track_length - 1)),
                continue;
            end
            
            %make a column of angle changes without non-motile cell angles as
            %defined by comparing the average dist for a given cell to the
            %user's entered value
            
            avg_dist = mean(del_dist(start_pos:end_pos));
            motile2 = [];
            
            if (avg_dist >= min_disp),
                motile2 = del_angle(start_pos:end_pos-1);
            else
                motile2 = NaN*ones(size(del_angle(start_pos:end_pos-1)));
            end
            
            %print first line of analysis output
            fprintf(analysis_fid, '%d\t%d\t%5.2f\t%5.2f\t%5.2f\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%5.2f\t%s\t%5.2f\n', cur_data(1,1), 1, cur_data(start_pos, 3:5), ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',cur_data(1,2), ' ',cur_data(1,6));  %modified 2/19/2011 to make MATLAB consider empty cells as such
            %print second line of analysis output (no angle)
            fprintf(analysis_fid, '%d\t%d\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%s\t%5.2f\t%s\t%s\t%5.2f\t%5.2f\t%5.2f\n', ...
                cur_data(1,1), 2, cur_data(start_pos+1,3:5), del_dist(start_pos), del_pos(start_pos,:), inst_vel_3D(start_pos), inst_vel(start_pos,:), del_disp_3D(start_pos), del_disp(start_pos,:), ' ', inst_vel_2D(start_pos), ' ', ' ', cur_data(start_pos+1,2), del_time(start_pos), cur_data(start_pos+1, 6)); %modified 1/27/2011 to make MATLAB consider empty cells as such, on 2/19/11 it also reports cell length
            %print the rest of the analysis file
            for t = (start_pos + 1):end_pos,
                fprintf(analysis_fid, '%d\t%d\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\n', ...
                    cur_data(1,1), t+1, cur_data(t+1,3:5), del_dist(t), del_pos(t,:), inst_vel_3D(t), inst_vel(t,:), del_disp_3D(t), del_disp(t,:), del_angle(t-1), inst_vel_2D(t), motile(t-1), motile2(t-start_pos), cur_data(t+1,2), del_time(t), cur_data (t+1,6));
            end    
            fprintf(analysis_fid,'---\t\n');    
            
            xyz_del_disp = cur_data(end_pos+1,3:5) - cur_data(start_pos,3:5);
            ttl_del_disp = sqrt(sum(xyz_del_disp.^2));
            xyz_del_dist = sum(abs(del_pos(start_pos:end_pos,:)),1);
            ttl_del_dist = sum(del_dist(start_pos:end_pos));
            ttl_del_dist_2D = sum(del_dist_2D(start_pos:end_pos));  % line added 2/9/2011
            CI = ttl_del_disp/ttl_del_dist;
            xyz_avg_vel = mean(inst_vel(start_pos:end_pos),1); 
            ttl_avg_vel = mean (inst_vel_3D(start_pos:end_pos));
            avg_del_angle = my_nanmean(del_angle(start_pos:end_pos-1));
            avg_2d_vel = mean(inst_vel_2D(start_pos:end_pos));
	        duration = sum(del_time(start_pos:end_pos)*cyc_time/60);  % LINE ADDED 1/16/2011
            mean_vel_2D = ttl_del_dist_2D/duration;  %line added 2/9/2011
	        mean_vel_3D = ttl_del_dist/duration;  % LINE ADDED 1/17/2011
	        arr_coeff = sum(nonmotile_arr_coeff(start_pos:end_pos)) / sum(del_time(start_pos:end_pos)); % LINE ADDED 1/17/2011
            
            
            %print out summary line
            fprintf(summary_fid,'%d\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f', ...
                cur_data(1,1),ttl_del_dist, xyz_del_dist(1:3), mean_vel_3D, xyz_avg_vel(1:3), ttl_del_disp, xyz_del_disp(1:3), avg_del_angle, CI, avg_2d_vel, duration, arr_coeff, mean(cur_data(:,6)));
            if length(exceeded)>2,
                fprintf(summary_fid,'\t*');
            end
            fprintf(summary_fid,'\n');
            
            % calculate and print out msd
            for step = 1:(end_pos-start_pos+1),
                dist = [];
                for cur_pos = start_pos:(end_pos+1-step),                
                    dist(cur_pos-start_pos+1) = sqrt(sum((cur_data(cur_pos+step, 3:5) - cur_data(cur_pos, 3:5)).^2));
                end
                fprintf(msd_fid, '%d\t%d\t%5.2f\t%5.2f\n', cur_data(1,1), sum(del_time(1:step))*cyc_time, mean(dist), mean(dist)^2);  
            end
            fprintf(msd_fid, '-----\t\n');
        end
        waitbar(0.25 + 0.75*(i/length(all_data)));
    end
    
    varargout{1} = all_data;
    
    fclose(analysis_fid);
    fclose(summary_fid);
    fclose(msd_fid);
    close(wh);
% catch
%     fprintf('An error occurred during tracker file processing.\nLast Error: %s\n', lasterr);
%     fclose(analysis_fid);
%     fclose(summary_fid);
%     fclose(msd_fid);
%     close(wh);
% end

